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THE  ESTIMATION  OF  DELAY  GRADIENTS  FOR  PURPOSES 
OF  ROUTING  IN  DATA-CQMMUNICATION  NETWORKS* 
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Adrian  Segal*** 


ABSTRACT 

\) 

The  derivatives  with  respect  to  message  flow  of  the  total  delay 
accumulated  per  unit  time  on  each  link  in  a  Data-Communi cation  network 
have  been  shown  to  be  fundamental  quantities  in  the  solution  of  the 
routing  problem.  Casting  the  problem  of  estimating  these  delay  gradients 
in  a  queueing  theory  framework,  and  making  no  statistical  assumptions 
other  them  stationarity,  we  propose  three  algorithms  that  process  the 
record  of  arrivals  and  departures  of  a  single-server  queue  to  derive  an 
estimate  for  the  derivative,  with  respect  to  arrival  rate,  of  the  total 
delay  accumulated  per  unit  time.  Through  simulation  and  analysis  we 
show  that  all  three  algorithms  are  asymptotically  unbiased  and  efficient 
for  M/D/1  queues.  By  simulation  of  other  queues  we  investigate  the  rela¬ 
tive  robustness  of  the  three  procedures.  Finally,  through  examination 
of  the  storage  and  computational  requirements  we  identify  a  single  most 
promising  algorithm. 
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1.  INTRODUCTION 

In  several  recent  works  [l]-[6]  it  has  been  shown  that  routing  in  a 
data-communication  network,  where  the  goal  is  to  improve  and  minimize  net¬ 
work  average  delay,  should  be  based  on  knowledge  of  the  derivative  of 
delay  on  the  network  links  with  respect  to  the  link  flows.  To  be  precise, 
if  we  let  A.  be  the  average  message  flow  (in  messages/sec.)  on  a  link 
(i,k)  of  a  communication  network  and  assume  that  the  average  delay  faced 
by  all  messages  going  through  this  link  is  a  function  only  of  A^,  then 
we  can  write  for  the  total  delay/unit  time  on  the  link,  namely 

the  product  of  the  average  delay  per  message  and  the  average  flow  of 
messages/unit  time.  With  these  notations,  several  of  the  centralized 
and  distributed  routing  algorithms  introduced  in  the  above  mentioned 
references  employ  the  derivative  D^(A  =  dD^/dA , ^  as  a  basic  quantity. 
One  way  to  obtain  this  quantity  is  to  make  some  statistical  assumptions 
on  the  message  interarrivals  and  lengths  and  use  one  of  the  queueing 
formulae  to  obtain  an  explicit  relationships  between  and  A^.  For 
example,  for  an  M/M/1  queue  one  obtains  (A^  1  *  -  V 

and  now  by  measuring  the  flow  A  ,  the  bit  capacity  C  ,  and  1/U..  ,  the 

iX  IX  3-X 

average  message  length  in  bits,  one  can  find  (A.  •  Of  course,  dif¬ 
ferent  statistical  assumptions  will  give  different  answers  for  > 

and  in  fact  these  explicit  formulae  can  be  obtained  only  for  the  situation 
where  independence  and  identical  distribution  of  the  variables  involved 
is  assumed.  To  circumvent  some  of  these  problems  a  different  approach 
can  be  taken.  Rather  than  estimate  the  flow  and  substitute  into  sone 
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explicit  formula  for  Djjc  (X^) ,  we  can  instead  f: algorithms  to  directly 

estimate  D*  (X  ) ,  the  latter  being  independent  of  any  statistical  assump- 
iiC  iK 

tions.  The  purpose  of  this  paper  is  to  introduce  and  analyze  three  proce¬ 
dures  for  estimating  the  derivative  D* (X)  of  the  total  delay  D(X)  on  a 
link  as  a  function  of  the  message  flow  rate  X  on  this  link.  The  algorithms 
are  based  on  measurements  of  arrival-departure  records  of  the  messages  and 
on  processing  these  quantities  recursively  to  obtain  the  required  estimates. 
We  use  the  nomenclature  of  queueing  theory  by  referring  to  the  link  as  a 
queue  and  to  the  messages  as  customers. 

In  Section  2  we  introduce  and  motivate  the  three  procedures  and  pre= 
sent  their  realization  in  recursive  form.  Section  3  sketches  the  analysis 
of  the  statistical  properties  of  each  algorithm  as  applied  to  M/D/1  and 
M/G/l  queues.  In  Section  4  we  compare  and  contrast  the  computational  and 
storage  requirements  of  each  procedure.  Finally ,  Section  5  presents  the 
results  of  simulating  the  algorithms  for  M/D/1,  M/M/1,  D/M/1 ,  and  D/M/1 
queues,  and  identifies  a  single  most  promising  procedure. 


2.  THREE  ESTIMATION  ALGORITHMS 

We  now  proceed  to  propose  and  evaluate  three  algorithms  which  process 
the  record  of  a  single  server  queueing  system  to  estimate  the  derivative 
of  the  total  delay/unit  Mm-  with  respect  to  arrival  rate  X.  The  available 
record  consists  of  exact  knowlecge  of  arrival  times  of  customers  to  the 
queue  and  of  their  departure  times  after  service  is  completed.  Time  is 
segmented  into  alternative  (random)  intervals,  busy  periods,  during  which 


I 


-3- 


the  server  is  occupied ,  and  idle  periods,  when  the  server  is  free.  The 

observation  interval  which  is  used  to  form  our  estimate  consists  of  a 

number  of  busy  periods  and  the  intervening  idle  periods. 

A  simple  thought-experiment  motivates  all  three  estimation  algorithms. 

Consider  our  single-server  queueing  system  with  its  average  arrival  rate 

of  X  customers  per  unit  time.  For  a  given  observation  period  T  ,  if  we 

E 

nan  compute  the  total  system  S,  i. e . ,  the  sum  of  all  customer's  service 

and  waiting  times,  then  the  average  delay /unit  time  is  given  by  D  *  S/T  . 
Suppose  we  could  actually  alter  the  input  flow  by  some  <5X .  Then,  on  the 
basis  of  an  earlier  D,  by  computing  D*  for  the  next  observation  period, 
we  can  estimate  the  derivative  of  the  total  delay/unit  time  by  calculating 


fi' 


D*-D 

<5X 


(2.1) 


However,  in  any  actual  queueing  system  It  would  be  undesirable  to  change 
flows  just  for  measurement  purposes.  Even  if  we  could  .'.mplement  (2.1)  , 
the  independent  statistical  fluctuations  in  D  and  D*  would  probably  make 
it  a  very  poor  estimator.  Hence,  what  we  need  is  some  mathematical  for¬ 
malism  for  an  imaginary  increment  in  flow  6X ,  which  will  allow  us  to  com¬ 
pute  the  corresponding  change  in  delay  without  actually  perturbing  the 
arrival  rate. 

According  to  intuition ,  an  increase  in  arrival  rate  should  result  in 

additional  customers  entering  the  system.  An  extra  customer  arriving  in 

a  time  interval  T  with  probability  £  will  increase  the  effective  rate 
E 

by  6X  »  e/T  .  If  extra  arrivals  are  mutually  independent  events,  the 
£» 


probability  of  two  or  more  customers  will  be  of  second-order  in  £  and 
hence  of  second-order  in  6X  .  Therefore ,  only  the  effect  of  a  single 
extra  arrival  has  to  be  considered  explicitly.  We  also  assume  that  the 
arrival  time  of  the  extra  customer  is  uniformly  distributed  over  the  ob¬ 
servation  period  T^.  In  addition,  in  order  to  explicitly  compute  the 
change  in  total  system  time  over  the  observation  interval  due  to  an  extra 
arrival,  we  will  assume  that  the  additional  customer  has  some  known  deter¬ 
ministic  service  requirement .  These  assumptions  allow  us  to  compute  an 
expected  increase  in  system  time  conditioned  on  the  arrival  of  a  new  cus¬ 
tomer,  and  the  resulting  estimation  procedure  will  be  called  the  customer- 
addition  algorithm. 

In  a  second  algorithm,  first  proposed  in  [1] ,  an  incremental  decrease 
in  the  effective  rate  X  is  simulated.  This  is  done  by  assuming  that  each 
customer  arriving  to  the  system  is  allowed  to  indeed  enter  the  queue  with 
probability  1  -  C,  and  is  eradicated  with  probability  e,  independently 
from  customer  to  customer.  In  this  way  we  simulate  an  arrival  process 
with  rate  X  (1  -  e)  ,  when  E  ■  -  5XA*  If  we  let  M'  be  the  total  number  of 
customers  in  a  period  T  ,  then  £  is  given  by 

L 


e  -  ^  6X  . 


(2.2) 


Again,  the  probability  of  removal  of  two  or  more  customers  from  the  same 

period  T  is  second-order  in  <5X  and  hence  the  reduction  of  total  system 
& 

time  that  has  to  be  considered  explicitly  is  due  to  removal  of 
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three  algorithms  will  occur  at  the  end  of  busy  periods,  and  assumes  the 
form 


O'  (k)  =  (k-1)  +  A^.  (2.4) 

Here,  D' (k)  denotes  the  delay  derivative  estimate  based  on  k  busy  periods 

A 

of  observation.  The  factor  ^  renormalizes  to  correspond  to  a  term 

in  the  expression  for  the  k-busy  period  estimate.  Finally,  A  denotes  the 

k 

contribution  to  the  k-busy  period  estimate  derived  from  the  current  busy 
interval. 

In  the  customer-addition  algorithm  we  simulate  an  increase  of  5X  in 
the  arrived,  rate.  The  following  assumptions  will  be  made: 

1)  The  probability  of  an  extra  arrival  in  the  interval  T  is 

£ 

6X*T  . 

E 

2)  Each  extra  arrival  is  independent  of  all  other  arrivals . 

3)  The  extra  arrival  is  uniformly  distributed  over  the  interval  T^,. 

4)  The  service  requirement  of  the  extra  customer  is  known;  we  denote 
it  by  x. 

We  must  note  at  the  start  that  the  form  of  the  customer-addition  algorithm 
derived  here  is  only  meaningful  for  queues  where  all  customers  have  identical 
service  requirements .  However,  the  procedure  can  be  practically  extended 
to  the  case  where  the  service  time  distribution  allows  a  finite  set  of 
discrete  service  requirements.  For  simplicity,  we  treat  here  the  case  of 
a  single  allowable  service  requirement. 
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Let  and  1^  denote  the  busy  and  idle  periods,  respectively. 

Let  denote  the  increase  in  system  time  over  N  busy  periods 

associated  with  the  arrived,  of  an  extra  customer.  We  let  6§  denote  the 

(N) 

expected  increase  in  system  time  associated  with  an  increase  in  arrival 
rate  6X  and  conditioned  on  the  record  of  arrival. s  and  departures .  By  con- 


ditionining  on  the  random  arrival  time  t  being  in  each  T  and  I  we  can 

•  Jv 

compute  5§  .  to  first  order  in  SX  by 


(N) 

6s 


(N) 


N  Tk 

{  l  E(lSs(N)lt  s  Tk'  lueueing  record)  — 
k=l  E 


N  Ik 

+  ^  E(5S(N)lt  6  Ik'  queueing  record)  —  }  T^-SX  {25) 

k=l  E 


the  factor  (T _,6X)  outside  the  brackets  is  the  probability  of  an  extra  ar- 
rival.  The  factors  T^/T^  and  I^/T^,  represent  the  probabilities  of  t  being 
in  the  k-th  busy  and  idle  periods,  respectively.  This  is  a  consequence 
of  the  assumption  that  t  is  uniformly  distributed  over  Tg.  Since  we  are 
interested  in  the  derivative  of  the  total  delay/unit  time  with  respect  to 
flow  rate,  our  estimator  is  given  by 


(N)  te 


(2.6) 


We  focus  next  on  the  calculation  of  E(6s  .  |t  S  T  ,  queueing  record) 

(N )  JC 

and  E(6s^^|t  8  1^,  queueing  record).  These  expected  increments  in  system 
time  are  composed  of  the  average  effect  of  the  extra  arrival  on  existing 
customers  plus  the  average  system  time  of  the  additional  customer.  In 
considering  additional  arrivals  in  a  busy  period,  we  can  distinguish 
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be tween  the  effects  on  customers  in  that  busy  period  and  interactions  with 
succeeding  busy  periods.  First  we  examine  the  part  of  E(6s^j  |t  6  T^, 
queueing  record) ,  call  it  As that  comes  from  considering  the  k-th  busy 
period  in  isolation.  To  facillitate  discussion,  the  following  notation  is 
defined: 


a  arrival  time  of  i-th  customer  in  busy  period  relative  to 
start  of  busy  interval. 

=  service  requirement  of  the  i-th  customer 
*  system  time  of  the  i-th  customer 
M  =  number  of  customers  in  the  busy  period 
T  =  duration  of  busy  period. 


The  system  time  of  a  given  customer  is  equal  to  his  service  plus  waiting 
time.  The  waiting  time  is  equal  to  the  sum  of  the  service  requirements 
of  those  who  enter  the  busy  period  before  him  minus  his  arrival  time. 
Hence,  for  an  arbitrary  queue,  the  system  time  of  the  i-th  customer  is 
given  by 


s. 

1 


X,  + 

i 


i-1 


T. ) 
a 


(2.7) 


Now  we  note  that  an  additional  arrival  at  time  t  in  the  interval 

i 

(T  (T,  I  will  have  to  wait  (  £  x,  -  t)  until  being  served  and  will  cause 
1  1+1  £=1  * 

the  (M  -  i)  following  customers  to  wait  an  additional  time  x.  Assuming 
that  all  the  service  requirements  are  the  same  (x^  =  x) ,  an  arrival  at 
time  t  results  in  an  additional  total  waiting  time  of  (M  -  i)x  +  ix  -  t  or 
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Mx  -  t.  Now  letting  t  be  a  uniformly  distributed  random  variable  on  0  to 
T  or  equivalently  on  0  to  Mx,  we  note  that  t  =  T/2  or  t  =  Mx/2.  Hence, 
the  expected  increment  in  system  time  As^  is  given  as 


Ask  -  x  +  4»k  , 


(2.S) 


where 


Aw^  =  Mx/2 


(2.9) 


We  complete  the  calculation  of  E(6s  t  §  T ,  queueing  record)  by 

IN]  JC 

considering  the  additional  system  time  that  may  result  from  one  busy  period 
overlapping  onto  another.  No  matter  where  an  additional  customer  arrives 
in  the  k-th  busy  period,  that  period  will  be  extended  by  the  service  time 
x.  The  value  of  x  relative  to  the  following  idle  period  durations  will 
determine  the  number  of  succeeding  busy  periods  that  will  be  affected  by 
an  arrival  in  T^.  If  x  <  no  following  busy  periods  suffer  additional 
delay.  If  I  <  x  <  I  +  I  only  the  (k-t-l)-st  busy  period  is  affected. 

K  ~ ~  JC  JCi"  X 

The  exact  effect  on  a  given  busy  period  j  in  the  future ,  depends  on  how 
much  an  arrival  in  T^  causes  the  (j-l)-st  busy  period  to  overlap  into  the 
j-th  busy  period.  For  example,  if  x  >  1^,  then  each  customer  in  T^+^  will 
suffer  an  additional  delay  (x  -  1^) .  Letting  denote  the  number  of  cus¬ 
tomers  served  in  the  k-th  busy  period,  the  preceding  reasoning  leads  to 

the  following  rule  for  computing  E(Ss  |t  6  T  ,  queueing  record) : 

C.N]  K 
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(2.10) 


To  complete  our  description  of  the  customer-addition  algorithm,  we 

just  now  evaluate  the  average  increase  in  system  time  E(5s  jt  S  I  , 

\  N  / 

[ueueing  record)  associated  with  arrivals  in  idle  periods.  The  effect  of 
m  arrival  in  1^  on  the  (j+k)-th  busy  period  again  depends  on  how  much  the 
[j+k-l)-st  busy  period  slides  onto  the  (j+k)-th  busy  interval,  and  thus 
just  be  averaged  over  all  times  of  arrival  t  of  the  additional  customer  in 
die  k-th  idle  period.  Let  t'  denote  the  negative  of  the  time  between  the 


irrival  of  the  additional  customer  and  the  start  of  the  (k+l)-st  busy  inter- 


ral.  Then,  assuming  an  additional  customer  arrives  in  1^,  t'  is  a  random 

rariable  uniformly  distributed  on  the  interval  [-1^,0].  Let  us  define 

die  quantity  O  .by: 

k,  j 


m®! 


k+m 


1-1 

3*1 


(2.11) 
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Ti.en,  an  additional  arrival  in  I  will  cause  the  (j+k-1)  -th  busy  interval 

to  overlap  onto  the  ( j+k) -th  busy  period  by  an  amount  equal  to 

max {x  +  t'  -  a  . ,  0}.  Hence,  each  customer  in  the  (j+k) -th  busy  period 
k,3 

suffers  an  additional  waiting  time  of  max{x  +  t'  -  a . ,  0}.  On  the  basis 

k,  j 

of  this  reaisoning,  we  may  compute  E{6s  |t  6  I  ,  queueing  record}  by 

(N )  K 


4+1 


^  , 

E{6s  .  |tei.,  queueing  record}  =»  x  +  £  M.  .  J  max{x+t'-a  .  ,0}r — 
(NJ  k  j-i  K  V=-I.  k,D  * 


for 


4  4+1 

A,  Vj  -  X  -  A,  ^k+j 


j-1 


j=l 


(2. 12) 


(2.13) 


The  integrals  appearing  in  relation  (2.12)  are  computed  as 

I  \+^x — -  "  J  (-lu  >  a,,  ,-x) 

Vi,/.  -  »  * 

^  ^  (x-a  .)• 

k  k'3 


r. 


s+j 

Vi 


k#  j  ‘ 
2 


k,j 


(a,  ,-x  >  -i  ) 

k,j  k 


(2.14) 


Employing  relations  (2.8),  (2.10),  and  (2.12)  in  (2.5)  and  (2.6)  we 
can  conceive  of  a  processor  which  updates  an  estimate  for  the  delay  gradi¬ 
ent  at  the  end  of  each  busy  period.  Let  4^  denote  the  time  from  the  start 
of  the  observation  period  to  the  end  of  the  k-th  busy  interval.  Let 
denote  the  incremental  expected  waiting  time  suffered  by  the  M,  customers 
in  the  most  recent  busy  period  due  to  an  additional  arrival  during  the 
entire  current  queueing  record  plus  the  expected  whiting  time  of  an  extra 


customer  arriving  in  T^. 


Then,  the  estimate  updating  at  the  end  of  the 


k-th  busy  interval  assumes  the  form 

.-»>  *  \ 


3s=i  & 

\  *-11 


(2.15) 


the  (D'  -x)  terms  appear  since  we  can  identify  (D '  -x)  as  am  estimate 

(il  U) 

of  the  incremental  increase  in  waiting  time  per  unit  time. 

Now,  A^  can  be  expressed  as 


Aw,  +  A 


\  "Ik 


(2.16) 


where  Aw^  denotes  the  expected  additional  waiting  time  suffered  by  members 

of  the  current  busy  period  due  to  an  additional  arrival  in  plus  the 

newcomers  expected  wad. ting  time,  while  A  and  A  ,  denote  the  expected 

Ik  BX 

waiting  time  suffered  by  customers  in  the  most  recent,  k-th,  busy  inter¬ 
val  due  to  an  arrival  in  either  past  idle  periods  or  busy  periods  respec¬ 
tively.  Since  x  is  finite,  we  need  only  look  back  a  finite  number  of  busy 
and  idle  periods  to  compute  A^  and  Let  n^  denote  the  index  of  the 

first  idle  period  at  which  a  new  arrival  with  service  requirement  x  can 
affect  the  customers  of  the  current,  k-th,  busy  interval.  We  determine 
n^  by  requiring  that 


k-1  k-1 

l  li<  x<  l  h 

^V1  ^k 


(2.17) 


Then,  we  may  employ  relations  (2.10)  and  (2.12)  to  compute  A  amd  A  by 

BX  XX 
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k-1  k-1  T 

^Bk  *  ^  “k1*  '  .  £  V  ^  ' 

Jjptr  +1  J  =  z  J  K 


(2.18) 


and 


Ik 


?k. ./. 


*™k  1 


max{x+t'-a„. 


0> 


dt’  )  X4 

■  -rr— 


t'— I, 


Jt,k— flJL  | 


T 

k 


(2.19) 


We  note  that  at  the  end  of  each  busy  interval,  n^_  must  be  updated  so  that 
relation  (2.17)  is  satisfied. 

In  the  customer-removal  algorithm  we  simulate  a  decrement  in  arrival 
rate  6A  by  removing  customers  from  the  queue  with  probability  e  (given  by 
relation  (2.2))  and  computing  the  resulting  decrement  in  total  system  time. 
Here,  the  situation  is  somewhat  simpler  than  in  the  customer-arrival  algo¬ 
rithm,  since  removal  of  a  customer  affects  only  the  customers  belonging  to 
the  same  busy  period.  Letting  6s  ^  denote  the  change  in  the  system  time 
of  the  i-th  busy  period  due  to  the  removal  of  the  j-th  customer,  the  ex¬ 
pected  change  in  system  time  6s  in  a  N  busy  period  observation  interval 

(NJ 

due  to  a  decrement  in  flow  6A  and  conditioned  on  the  queueing  record 


is  formulated  as 


N  i 


Te6A 


E  (6s I  queueing  record)  =  V  7  -  6s  ,  (  ) 

(N)  i-1  j-1  3'1  M 


(2.20) 


Hence,  the  desired  delay  gradient  estimator  is  given  by 

M 


A  .  E  (,6s  |  queueing  record) 

n1  =  — - iiii -  . 

(N)  T£  6A 


1  N  i 

ffr  I  l  -  *«.  t  • 

i-1  j-1  3,1 


(2.21) 


We  compute  6s.  .by  working  with  more  microscopic  quantitites.  Let  C 

j  r  i  in  t  i 
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denote  the  amount  of  system  time  saved  for  the  m-ch  customer  in  the  i-th 
busy  period  by  the  removal  of  the  n-th  customer  in  that  busy  period.  Since 


the  removal  of  the  n-th  customer  can  have  no  effect  on  customers  that  pre¬ 
ceded  him,  Cn  .  *  0  for  m  =  l,...,n-l.  Hence,  5s.  .  can  be  computed  as 
m,i  j,r 


•  -  -  I  cj 

I'1  “'i 


(2.22) 


We  now  develop  a  systematic  procedure  for  calculating  the  ^ ' s . 

To  simplify  notation,  we  drop  the  i  denoting  the  index  of  the  busy  period. 
Letw^,  s^,  and  denote  the  waiting  time,  system  time  and  service  require¬ 
ment,  respectively  of  the  n-th  customer  in  the  busy  period.  Let  dn  and 
an  denote  the  corresponding  departure  and  arrival  time  of  the  n-th  customer. 

Since  the  system  time  the  n-th  customer  saves  by  its  own  removal  is  s  , 

n 

we  have  Cn  =*  s  .  In  considering  the  effect  of  removing  the  n-th  customer 
n  n 

on  the  (n+l)-st  customer,  either  a  new  busy  period  begins  with  the  (n+1) -st 

or  the  (n+1) -st  customer  remains  part  of  the  busy  period  formed  by  customers 

1  to  n-1.  The  condition  for  customer  n+1  beginning  a  new  busy  period  is 

that  the  arrival  time  a  of  the  (n+l)-st  customer  is  greater  than  the 

n+1 

departure  time  d  of  the  (n-l)-st  customer.  In  this  case,  customer  n+1 
n-l 

will  save  its  waiting  time  w  . , .  If  d  ,  >  a  , ,  removal  of  customer  n  does 

n+1  n-l  n+1 

not  cause  customer  n+1  to  start  a  new  busy  period,  and  hence  customer  (n+1) 


(2.23) 


Examination  of  the  definition  (2.23)  for  cn  ,  reveals  that  it  may  be  more 

n+l 


succinctly  stated  as 


Vl  * 


(2.24) 


Similar  reasoning  to  that  employed  in  calculating  cQ+i  applies  to  the 

computation  of  c11.  The  removal  of  customer  n  either  causes  customers  m  and 

m-1  to  be  in  the  same  busy  period,  or  breaks  the  busy  period  after  customer 

(m-1) .  The  removal  of  customer  n  causes  customer  m-1  to  3ave  system  time 

Cn  Hence,  customer  m-1  departs  at  an  earlier  time  d  ,  -  c°  ,  .  If 

m-1  n-1  m-1 

this  new  departure  time  for  customer  m-1  is  greater  than  the  arrival  time 


a  of  customer  m,  customers  m  and  m-1  remain  in  the  same  busy  interval 
m 

and  customer  m  saves  a  system  time  c”  , .  However,  if  a  >  d  -  cn  , 

m— 1  m  m-1  m— 1 

•customer  m  begins  a  new  busy  period  and  saves  its  waiting  time  w  .  These 

m 

relationships  are  summarized  in  the  following  rule  for  computing  c11. 

m 


|  Vi 


for  d  ,  -  c  >  a 
m-1  m— 1  —  m 


'a  I  .  ,  _  .  n 

Id  -  a  ■  w  for  a  >  d  -  c 

mm  m  m— 1  m— 1 


(2.25) 


This  rule  may  be  expressed  more  compactly  as 


c11  *  min{cn  ,  ,w  }  . 
m  m— 1  m 


(2.26) 


Hence,  the  algorithm  for  computing  the  c  ’s  may  be  summarized  as 

m 


In  our  third  algorithm,  we  simulate  am  increase  in  rate  6X  by  a  lineau: 
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contraction  in  time  scale.  Assume  that  t.  .  denotes  the  arrival  time  of 

J  ** 

the  j-th  customer  in  the  i-th  busy  period  relative  to  the  beginning  of  the 
observation  interval.  We  define  a  new  set  of  shifted  arrival  times  by 


t*  ■  (l  - 

j,i  U  X  j,i 


(2.32) 


Since  iSX  represents  an  infinitesimal  change  in  rate,  we  can  choose  it  suf¬ 
ficiently  small  so  none  of  the  busy  periods  are  shifted  onto  other  busy 
periods  by  the  time  contraction. 

We  now  consider  the  increment  in  system  time  that  comes  from  each  cus¬ 
tomer  arriving  a  little  earlier  and  hence,  waiting  a  little  longer.  The 
waiting  cime  of  the  j-th  customer  in  the  i-th  busy  period  is  defined  by 


w 


- 3? « 

J-1  nil  ^ 


(2.33) 


where  x.  denotes  the  service  requirement  of  the  Z-th  customer  in  the 

i-th  busy  period.  If  we  substitute  the  shifted  arrival  times  given  by 

(2.32)  into  (2.33),  we  can  relate  the  new  waiting  times  w*  .  to  w .  by 

j/i 


■  w 


j»i  j.i 


♦&t 


j»i 


- 


(2.34) 


We  define  r .  .  as  the  arrival  time  of  the  j-th  customer  in  the  i-th  busy 
j 

period  relative  to  the  start  of  that  busy  period.  Hence,  the  additional 
system  time  over  N  busy  periods  that  follows  from  (2.34)  can  be  expressed 
as 


f  l  \ ' 

A  i-1  1 


(2.35) 


where 


"i 

l  T*,i 


M.  -1 


M  >  1 


(2.36) 


A  second  contribution  to  the  increment  in  system  time  follows  from 

t-he  fact  that  our  contraction  procedure  shifts  the  right  edge  of  the 

5A 

observation  interval,  leaving  a  gap  of  duration  - — T  in  which  additional 

A  E 

customers  could  arrive .  Let  D  (X)  denote  the  average  system  time  per  cus- 

c 

tomer.  Then,  XD  (X)  represents  the  average  total  delay  accumulated  per 
c 

unit  time.  Hence,  the  average  increment  in  system  time  associated  with 

5A 

an  interval  of  length  would  be  given  as 

A  E 


(t~t_)  (Ad  (X) )  -  6XT  d  (X)  . 

A  E  C  EC 


(2.37) 


Since  both  X  and  D  (X)  are  unknowns ,  we  use  the  fact  that  asymptotically, 
c 


TE  ' 


(2.38) 


m  M- 
N  i 


Dc(X)  - 


I  Iv 

i-1  j-1 


N 

l  Mi 

i-1 


(2.39) 


l 
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where  S.  .  denotes  the  system  time  of  the  j-th  customer  in  the  i-th  busy 
3  »x 

period.  Combining  the  contributions  to  the  incremental  system  time  given 

by  (2.35)  and  (2.37)  with  relations  (2.38)-(2.39)  we  obtain  the  following 

N-busy  period  delay  gradient  estimator: 

M. 


,  K1  k  8  J  *  zj 


IN) 


M 


(2.40) 


(N) 


where  Z.  is  defined  by  relation  (2.36)  and  by  (2.29) 

i  (N) 

may  be  written  in  recursive  form  as 

M. 


D’ 


(k) 


-ttiiL  6'  ♦  A 

M(k)  (k-1)  \ 


Relation  (2.40) 


(2.41) 


where 


4  r 

\  “  I  Sl,*  *  \ 


(2.42) 


1-1 


Alternatively,  by  employing  the  definition  of  S?  ^  and  Z^,  we  may  express 

\  “ 

^  (Mk-4  +  1Ul,k  * 


(2.43) 


3.  STATISTICAL  ANALYSIS  OF  THE  ALGORITHMS  FOR  M/D/1  AND  M/G/l  QUEUES 

Having  motivated  and  defined  all  three  algorithms  for  estimating  the 
delay  gradient,  we  now  proceed  to  analyze  their  asymptotic  behavior  as 
the  observation  interval  and  hence  the  number  of  busy  periods  N  included 


i 
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become  large,  for  the  special  cases  of  some  simple  queues.  We  sketch  the 
analysis  of  the  asymptotic  bias  for  all  three  algorithms  in  the  case  of 
an  M/D/1  queue.  In  addition,  for  an  M/D/1  queue  we  derive  Cramer-Rao 
lower  bounds  for  the  r.m.s.  error- of  any  unbiased  estimator  using  the  same 
information  as  each  of  the  three  algorithms.  Finally,  we  will  discuss 
briefly  the  analysis  of  the  asymptotic  bias  of  the  customer-removal  and 
time-contraction  algorithms  in  the  case  of  M/G/l  queues. 

For  a  general  queue,  the  asymptotic  bias  b  associated  with  each  algo¬ 
rithm  is  formulated  as 


where  D(X)  is  the  average  total  delay  accumulated  per  unit  time.  If  Dc(A) 


denotes  the  average  delay  per  customer,  then  D(A)  is  expressed  as  ADc(X), 
and  hence  fr-  may  be  written  as 


3D 

ax 


D  (X)  +  X 
c 


3D  (X) 
c 

3X 


(3.2) 


Now  the  average  system  time  per  customer  Dc(A)  is  expressible  as 

D  (X)  »  x  +  w(A)  , 
c 

where  x  denotes  the  average  service  time  and  w(A)  denotes  the  average 
waiting  time.  Hence,  |^-  can  be  reformulated  as  a  function  of  the  average 
service  requirement  and  waiting  time  as 
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3D  ,3w(X) 

»mx*  W(X1  *  xi — 


(3.4) 


We  can  evaluate  the  above  expression  for  all  queues  for  which  ah  explicit 
form  of  the  waiting  time  distribution  is  available.  For  an  M/G/l  queue, 
w(X)  is  given  in  [7]  as 


P(1  +  o 
D 

2(1  -  p) 


where 


p  *  Xx  , 


(3.5) 


(3.6) 


is  the  utilization  factor  and 


(3.7) 


where  a,  denotes  the  variance  of  the  service  time  distribution, 
b 

We  now  briefly  describe  the  steps  employed  to  show  the  asymptotic  un¬ 
biasedness  of  the  customer-addition  algorithm  in  the  case  of  an  M/D/1 
queue.  For  the  full  details  we  refer  the  reader  to  [8].  From  relation  (3.1) 
we  must  evaluate 


lim  E  6'  „  =  lim  E{^ - rr^->  , 


(3.8) 


where  <5s is  defined  by  relation  (2.5),  in  order  to  compute  the  asymptotic 
(N) 

bias.  Letting  M  denote  the  average  number  of  customers  per  busy  period,  by 
the  law  of  large  numbers  we  have 
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11m  T 

E 

N*» 


lim{^}  , 

r> 


C3.9) 


and  interchanging  the  limit  and  expectation  operations  in  (3.8)  we  need 
to  compute 


A  lim  {—  }  .* 

t.'-**1  NM  °A 


(3.10) 


Now  define  the  two  aggregate  vectors 


M 


(N) 


(3.11) 


and 


X(N)  =  (I1'*-,'IN-1)  ' 


(3.12) 


where  the  1 s  are  the  number  of  customers  served  in  the  i-th  busy  interval 


and  the  1^ 1 s  are  the  i-th  idle  period  durations .  Then  for  an  M/D/1 


queue,  the  queueing  record  is  completely  specified  by  M  and  I 


break  up  the  calculation  of  E{- 


<5S 


(N) 


(N)  * 


We 


(N) 


-tt — }  by  first  conditioning  on  M  , 
0a  (N) 


averaging  over  ,  and  performing  a  final  averaging  over  M ,  . 

(N;  (N) 

<5* , 


Writing 


out  E{- 


(N) 


<S\ 


}  from  relation  (2.5)  as 


-  X  e1e!6s<“'1  tsT,<’  “w  *  X  E{E,ss(»)|tsv,N)'i(») 


(N)  (N) 

(3.13) 


We  note  that  there  are  two  types  of  terms  that  must  be  analyzed,  corresponding 
to  the  busy  and  idle  period  contributions  to  the  incremental  delay.  To 
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V*1  ^ 

perform  the  Initial  averaging  of  — — —  over  I  we  use  our  expressions  for 

E(5S(N)  lt6V  *(N)'  *(N)}  E(5S(nJ  tSIk'  “(N)'  ^(N)5  given  by  relati°ns 

(2.10),  (2.12)  and  the  fact  that  the  I^'s  are  independent  exponentially 

distributed  random  variables  with  probability  densities 

PI  (T)  =  XeXT  °-i(T)  *  (3.14) 

j 


After  considerable  manipulations,  which  are  detailed  in  [8],  we  obtain  the 
5s , 


mean  of 


(N) 


conditioned  on  M  ^  as 


<Ss 

E{— M-l- 


N 


2  122  N-l  N-j 


«  '”<»>*  ’  *  2V  1  *  £ 


where 


N-l 


+  Ji  ‘*  +  K+i<x  +  x<e"AX  -  “>> 


1  -Ax 


N-2  N-k 

+  y  y  m,  .b. 

k=l  j=>2  k+*  3 
N-2  N-k 

+  y  y 

,  L.  .*•_  M  .c. 
k=l  j=2  k+j  j 


(3.15) 


k.-Ax  r  (Ax)* 


,  K.  “Aa  r 

ak  =  x(x  -  r)e  J 


X  4=»k  11 


+  x  e 


2-Ax  (Ax) 


k-l 


(k-1)  !  ' 


(3.16) 


b.  =.  zXx{$  y  _ 

3  (£+2) 1 


tt*)1+2_j_  j  , 


(3.17) 


and 


(Ax)  ^  -Ax  x 

c  *  -  ■  ■  •  e  — 

j  (j+l)l  A  ‘ 


(3.18) 


We  now  note  that  the  M^ ' s  are  independent,  identically  distributed  random 
variables  and  that  for  an  M/D/1  queue  the  first  and  second  moments  M  and 
M*  are  given  by 


(3.19) 


By  using  (3.15)  and  (3.20)  we  can  compute  the  unconditional  expectation  of 
<$S 

(N) 

— .  Then  by  employing  (3. 15) -  (3 . 18)  we  are  able  to  evaluate  the  limit 


lim  — 

N-xx>  vjm 


(3.22) 


The  above  result  (3.22)  is  exactly  the  quantity  obtained  by  specializing 
the  expression  (3.5)  for  w(A)  to  the  case  of  an  M/D/1  queue  and  substituting 
into  our  general  expression  (3.4)  for  the  delay  gradient  -  .  Hence,  the 

customer- addition  algorithm  is  asymptotically  unbiased  for  an  M/D/1  queue. 

Since  the  calculation  of  the  exact  variance  associated  with  the  customer- 
addition  algorithm  is  too  cumbersome,  we  derive  a  Cramer-Rao  bound.  If  we 
have  an  observation  vector  R,  a  parameter  a  we  want  to  estimate  and  a 
density  for  R  parameterized  by  a  -  p(R:a) ,  the  Cramer-Rao  lower  bound  for 
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the  vatriance  of  any  unbiased  estimator  &(R)  of  a  is  stated  as  follows: 


Var  (&(R)  -  a)  > 


C[-3  In  p(R:g)j 
3a2 


(3.23) 


In  our  case  we  identify  R  as  (M  ,  I  ,„.  ) .  As  a  product  of  the  "memoryless" 

(N)  (N) 

property  associated  with  the  exponential  interarrival  time  distribution  we 


can  show  that  the  1 s  and  I j 1 s  are  all  mutually  independent  so  that 


P(S(N)'i:(N))  13  9±Ven  aS 


N 


1=1  j-1 


(3.24) 


where  for  an  M/D/1  queue  Pr{M.  =  m. }  is  given  in  [7]  as 

m,  -l 

(m.  p)  m_.  p 

Pr{M.  =  m. }  = 


m.  1 
i 


-  i 

e 


(3.25) 


Hence,  pretending  that  p  is  the  parameter  of  interest,  we  can  express 


P(M(N)'  I(H):P)  aS 


p  (m  ,  /  i  , :  p) 

*  (N)  (N)  H 


N  m. 

n  — 

i=l 


m.  -1 
i 


m.  1 


N  \  N  N-l 

l  VN  1  -  l  V  N-l  -fi.  I  1  ■ 

1=1  ’ei=1  (f) 


(3.26) 


But  now  letting  y  denote  for  an  M/D/1  queue,  where  is  defined  in  rela¬ 
tion  (3.4),  we  cam  express  p  as  a  function  of  y  as  follows: 
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Hence,  we  can  form  p  (M  » I  .  :  y)  and  can  evaluate  the  bound  in  relation 

(NJ  (NJ 

(3.23)  as 


Var  (y  -  <f)  > 


2  2 
JL_2 _ 


(3.28) 


(N-l+p)  (1-p) 


Next  we  investigate  the  asymptotic  properties  of  the  customer  removal 
algorithm  in  the  case  of  an  M/G/l  queue  by  first  interpreting  the  terms 
in  the  estimator.  For  a  given  busy  period,  the  inner  two  summations  in 
relation  (2.28)  may  be  grouped  into  two  terms  representing  the  sum  of  all 
the  service  times  of  customers  in  that  busy  period  and  the  cumulative 
service  mp  saved  by  by  all  other  customers  due  to  the  removal  of  each  customer 
separately.  Hence,  the  N  busy  period  delay  gradient  estimator  may  be  ex¬ 


pressed  as 


M. 

i 


I  I  s  Ip 

i-1  j-1  1,1  .  i=l  1 


M 


(N) 


M 


(N) 


(3.29) 


where 


Pi 


M.-l  M. 

"I  l  «“  i 

n-l  m=*n+l  m'1 


M.  =  1 
r 


M.  >  1 
i 


(3.30) 


We  examine  the  asymptotic  behavior  of  the  mean  of  the  estimator  speci¬ 
fied  in  (3.29)  by  interchanging  the  expectation  and  limit  operation.  By 
appealing  to  the  law  of  large  numbers  the  limiting  form  of  the  estimator 
as  N  becomes  unbounded  is 


“ 
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lia  E  D1  »  E  lim  D'  -  D  (X)  +  R  , 
N-~>  W  N—  tN)  C  5 


(3.31) 


where  D  (X)  is  the  average  system  time  per  customer,  p  is  the  expectation 
c 

of  the  quantity  defined  in  (3.30) ,  and  M  is  the  expected  number  of  customers 

served  per  busy  period.  Employing  the  general  relation  (3.4)  which  expresses 

the  delay  gradient  as  equal  to  D  (X)  plus  ,  we  can  formulate  the 

c  o  A 

asymptotic  bias  of  the  customer- removal  algorithm  as 


b  *  lim  E  D*  -  §  =  2-  -  X^$-  . 
n-h»  <N)  m  3X 


(3.32) 


We  can-  break  up  the  calculation  of  p  by  conditioning  on  M  =  i  for  i  =  2, — 
and  hence  we  find 


P  =  l  Q(i)f.  , 

i=2 


(3.33) 


where  f ^  is  the  probability  that  i  customers  are  served  in  a  busy  period 
and  Q (i)  is  defined  as 


i-1  i 

Q (i)  -  l  l  E(c"|M=i)  . 

n=l  m=n+l 


(3.34) 


Now  for  an  M/G/l  queue  the  average  number  of  customers  served  per  busy 
period  is  given  in  [7]  by 


1-P  * 


(3.35) 


The  z-transform  for  the  probability  density  of  the  number  of  customers 


served  per  busy  period  is  described  in  [7]  by  the  following  functional  equation: 
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F(z)  =»  zB*(X  -  XF(.z)}, 


(3.36) 


where  B* (s)  is  the  one-sided  Laplace  transform  of  the  service  time  density 
and  F (z)  is  defined  by 


F  (z)  »  £  f  z11 


(3.37) 


n“l 


3w  (X ) 


We  may  now  employ  the  expression  for  w(X>  in  (3.5)  to  write  X1X 
x3w(X)  _  x  ’  -  'V4.9WV4.11  V4.1 


as 


3X 


(i-p,  f  a  *  A  l  pktl 

*  D  k=0  ^ 


(3.38) 


where  p  =  Xx  is  the  utilization  factor.  Hence,  for  an  M/G/l  queue  the 
asymptotic  bias  may  be  expressed  as 


b 


(1-P)  [  l  Q(n)f 

n=2 


x(l+c^) 

4 


l  (k+2) (k+l)pk+1]  . 
k=*0 


(3.39) 


When  p  may  be  expressed  as  a  power  series  in  p,  relation  (3.39)  expresses 
b  as  an  expansion  in  powers  of  p. 

We  now  comment  briefly  on  the  calculation  of  the  Q(i) 's  in  relation 
(3.34).  A  complete  discussion  of  this  problem  is  contained  in  [8].  To 
calculate  Q(i),  we  need  to  evaluate  conditional  expectations  of  the  form: 

E{c"+k|M-i}  ,  (3.40) 


for  n  »  2...M,  k«l..M-n.  From  relations  (2.27),  we  can  express  cn+k  as 


cn  ,  *  min{x  ,w  „,w  •  .w  } 

n+k  n  n+1  n+2  n+k 


(3.41) 


I 


i 
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Our  approach  is  to  derive  the  density  for  (x  ,w  ..w  )  conditioned  on 

n  n+l  n+k 

M  =»  i  and  then  compute  F  (x)  ,  the  distribution  function  for  cn  ,  given 

n  n+k 

„  cn+k |M*i 

H  a  i  as  1 


F  n  (T)  -  1  -  Pr{*n  >  T.  »n+1  >  T,  »n+k.  >  T> 
cn+k  |M=»i 


(3.42) 


Now  we  calculate  P  (x)  as 

n 

Cn+k |M=*i 


P  n  (X)  *  d_  F  n  (X) 
Cn+k |M**i  Cn+k |M»i 


(3.43) 


and  derive  the  conditional  mean  (3.40)  as 


E(c 


n+k 


lM“i)  *  n 


(x)  dx  . 


'n+k  M-i 


The  calculations  become  simpler  in  the  case ' of  an  M/D/1  queue  where 


(3.44) 


W  ■  zn,k>  ■ 


(3.45) 


with 


z  .  ■  min{w  , . . . ,w  .  } 
n,k  n+l  n+k 


(3.46) 


Here  we  calculate  the  joint  density  of  (w  , . .w  ,  )  conditional  on  M  =*  i 

n+l  n+k 

and  then  compute  F  (X)  ,  the  distribution  function  for  z  given  M  *  i 

Zn,k|M-i  n'k 


F  (x)  =*  1  -  Pr{w  >  X,..w  >  x} 

z  .  | „  .  n+l  n+k 

n,k 


as 


(3.47) 
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Given  the  distribution  function  for  z  .  conditioned  on  M  *  i  we  compute 

n,x 

the  density  for  z  .  conditioned  on  M  =  i  by 
n,k 


n,k  |M-i 


(X) 

n,k  |M=i 


(3.48) 


and  can  write  down  the  conditional  density  for  cn+J<_  explicitly  as 

P  (T)  *  (1  -  F  (x)  )5(t-  x)  +  P  (t)  U-  (x  -  T)  (3.49) 

Cn+k|M-i  Zn,k  |M=»i  n,k|M-i“ 

By  employing  the  representation  for  the  asymptotic  bias  (3.39) ,  the 

technique  for  computing  the  Q(i)'s  sketched  in  the  preceding  paragraphs, 

and  known  formulae  for  the  f . 's  we  were  able  to  show  that  for  an  M/D/1 

x 

queue,  the  asymptotic  bias  b  may  only  contain  terms  of  fourth  order  or 
higher  in  p.  We  suspect  strongly  however  that  the  coefficients  of  all 
powers  of  p  in  the  expansion  of  b  obtained  from  (3.39)  vanish  for  an 
M/D/1  queue,  although  we  were  not  able  to  prove  this  explicitly.  Analysis 
of  relation  (3.39)  for  the  case  of  an  M/M/1  queue  reveals  that  b  contains 
nonzero  terms  of  order  p  and  higher. 

Having  analyzed  the  asymptotic  bias  associated  with  the  customer- 
removal  procedure,  we  now  proceed  to  derive  a  Cramer-Rao  bound  for  the 
variance  of  any  unbiased  estimator  of  the  delay  gradient  that  works  with 
the  same  observations  as  the  customer-removal  algorithm  in  the  case  of  an 
M/D/1  queue.  Since  for  an  M/D/1  queue  the  service  requirement  is  a  deter¬ 
ministic  quantity  -  x,  the  observations  which  the  customer  removal  algo¬ 


rithm  employs  are  the  M^'s,  the  number  of  customers  served  in  the  i-th 

busy  period,  and  the  w,  .'s,  the  waiting  times  of  the  j  =  2,...,M  customers 

j  r1  1 
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in  the  i-th  busy  period.  Hence,  we  define  the  N  busy  period  observation 
vector  as 


(3.50) 


where 


w 


2»i* 


,  .w 
M 


iri 


(3.51) 


Since  waiting  times  and  the  number  of  customers  served  per  busy  period 

are  statistically  independent  from  one  busy  period  to  the  next,  the  joint 

density  P (Y  .  )  may  be  expressed  as 
(N) 

N 


P(Y 


II  P(Y.) 
i-1  ^ 


(3.52) 


where  we  can  decompose  P(Yj  as 


P  (Y , )  -  P  (M.  )  P  (w  ..w  |  M.  ) 

i  i  2 ,  i  M .  .  i 


(3.53) 


i»l 


Employing  the  expression  for  p(M, )  given  by  (3.25)  and  the  expression  for 

p(w  ..w  I M. )  derived  in  [8],  we  can  form  P(Y  :p)  as 
*  2,1  M.  .  1  i 

9  «  i 


1,1 


i=l 


(3.54) 


(  l  M  -  N) 
i-1  1 


where  0  5.  w2  i  —  x  ' 


(3.55) 


0  <  w,  ,  .  <  x  +  w,  . 
-  k+l,i  -  k,i 


2. •Mi«l 


and 


(3.56) 
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Now,  using  the  relation  between  the  utilization  factor  p  and  the  delay 
gradient  y,  for  an  M/D/1  queue,  given  by  relation  (3.27),  we  may  form 


p(Y^:y)  and  compute  the  Cramer-Rao  band  from  (3.23)  as 


,  ~  2  x  p 

Var  (y  *  Y)  >, - y — T  • 

N(l-p) 


(3.57) 


Finally,  we  now  determine  the  asymptotic  behavior  of  the  time- 
contraction  procedure.  From  relation  (2.40),  we  find  by  exchanging  limit 
anrf  expectation  operations  and  appealing  to  the  law  of  large  numbers  that 


lim 

N-*» 


(3.58) 


where  z  denotes  the  expected  value  of  z^  defined  by  relation  (2.36). 
Hence,  by  the  same  reasoning  as  for  the  customer- removal  algorithm,  the 
asymptotic  bias  b  is  given  as 


(3.59) 


Similarly,  we  may  compute  z  by  conditioning  on  M  =  i,  i  =  2..°°  as 


z  -  l  G(i)f.  ,  <3-60> 

i-2 


where  f.  denotes  the  probability  of  i  customers  being  served  in  the  i-th 
busy  period  and  fi(i)  is  defined  as 

M 

ft(i)  -  E(  l  T  |M  =  i) 
i=2 


(3.61) 
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If  we  let 


(3.62) 


we  show  in  [8]  that  it  is  possible  to  derive  the  joint  density  of  the 

Tj's  for  j  ■  2..M  conditioned  on  M  *  i  and  hence  express  the  desired 

conditional  expectation  (3.61)  as  a  single  integral.  Thus,  in  the  same 

manner  as  for  the  customer-removal  algorithm,  relation  (3.59)  becomes  a 

representation  in  powers  of  p  for  the  asymptotic  bias  in  the  case  of  M/G/l 

9w  (X) 

queues.  By  employing  relation  (3.38)  for  we  obtain 

2. 

<30  x  (  1  +  C  )  00 

b  »  (1  -  p)  (  l  (I(i)f. - —  l  (k+2)  (k+1)  pk+X)  (3.63) 

i-2  1  4  k-0 


By  employing  relation  (3.63),  we  showed  that  for  am  M/M/1  queue  the  asymptotic 
bias  for  the  time-contraction  algorithm  contains  nonzero  terms  of  first 
order  in  p  amd  higher.  By  analyzing  the  form  of  the  time-contraction 
estimator  obtained  by  using  relation  (2.43)  to  define  we  are  able  to 
3how  asymptotic  unbiasedness  of  the  algorithm  for  an  M/D/1  queue.  From 
relation  (2.43),  we  may  alternatively  express  the  N  busy  period  estimator 


as 


N 


(3.64) 


where 


3 

X  (fl. 

3 


l+1)xU 


(.3.65) 


Hence  for  an  M/D/1  queue,  D'  would  be  given  by 


N  , 

Z  jK  +  V* 

j-1  3  -1 


(3.66) 


Now,  by  exchanging  limit  and  expectation  operations  and  appealing  to  the 
law  of  large  numbers  we  can  show  that 


~  1  ,m2  i 

lim  d;  -  r*{—  +  1} 

W  2  S 


(3.67) 


Employing  earlier  (3.19)  -  (3.20)  expressions  for  M  and  M,  we  can  verify 

3D 

that  the  limit  in  (3.67)  is  exactly  for  an  M/D/1  queue. 

We  now  conclude  our  analysis  of  the  time-contraction  algorithm  by 

deriving  a  Cramer-Rao  band  on  the  variance  of  any  unbiased  estimator 

of  the  delay  gradient  that  employs  the  same  observations ,  in  the  case  of 

an  M/D/1  queue.  Relation  (3.66)  indicates  that  for  an  M/D/1  queue  the 

observation  vector  consists  of  (M,  , . . . ,M„)  -  the  numbers  of  customers 

1  N 

served  in  each  busy  period  respectively.  Since  the  number  of  customers 
served  is  independent  from  one  busy  period  to  another  and  for  an  M/D/1 
queue  the  probability  of  M  customers  in  a  busy  period  is  given  by  (3.25) , 


p  ( M , . . M  : P )  is  given  as  follows: 

Mi-1 

N  (M. P)  M.P 

-  n  — ±—  ;  1 

.  M.  ! 

1-1  i 


p()«  ,..M  iP) 


(3.68) 
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Now,  by  employing  relation  (3.27) ,  which  expresses  p  as  a  function  of 

the  delay  gradient  v,  we  may  form  p(M  .  .M  :v)  and  hence  obtain  the 

IN' 

following  lower  band  on  the  variance  of  any  unbiased  estimate  y : 


VarCy  -  <f) 


2 

>  x 


fl 

N(l-p)5 


(3.69) 


4.  COMPUTATIONAL  COMPLEXITY  AND  STORAGE  REQUIREMENT  ANALYSIS  FOR  THE 
THREE  ALGORITHMS 

Having  completed  our  examination  of  the  statistical  properties  of 
all  three  estimation  algorithms •  we  proceed  to  compare  and  contrast  the 
computational  complexity  and  storage  requirements  of  each  procedure. 

Since  the  estimate  updating  at  the  end  of  the  k-th  busy  period  assumes 
the  form  in  relation  (2.4  )  for  each  algorithm,  we  consider  the  storage 
and  computational  requirement  of  forming  A^. 

For  the  customer  addition  algorithm,  A^  is  defined  by  relation  (2.16)  . 
To  compute  A  and  A_  we  need  to  have  stored  (I  ...I  )  and  T  ...T  ). 

\  A  “k  k_1  Hfcfl  k 

Letting 

"b,*  *  *  "  \  ,  (4.1) 

we  need  2n  storage  locations  for  the  busy  and  idle  period  durations 
b,k 

during  which  an  additional  arrival  can  cause  customers  of  the  most  recent 

k-th  busy  period  to  suffer  extra  delay.  Since  n  is  a  random  variable, 

b ,  k 

in  designing  a  system  to  implement  the  customer- addition  procedure  we  would 


(.4.2) 


a 


need  a  rational  for  choosing  a  buffer  size  N  (e) ,  such  that 

b 


Pr^^  >  Nfa(grj  }  <  e 


for  some  specified  tolerable  buffer  overflow  probability  e*  For  M/G/l 

queues,  a  manner  for  choosing  N,(e)  is  discussed  in  [8].  The  remaining 

b 

component  of  A^  is  A*^  which  is  defined  by  relation  (2.9).  Thus,  the 

storage  requirement  for  forming  A,  consists  of  2n  +0(1)  storage 

k  b  ,k 

locations,  where  0(1)  denotes  some  constant.  If  we  consider  the  compu¬ 
tations  involved  in  forming  Aw  ,  Al^r  and  ABk  we  count  a  total  of 


7nb,h+0(1) 


3"b,k 


+  0(1) 


2n  k  +  0(1) 

b ,  *• 


and 


b,k 


additions , 
multiplications , 
divisions , 
comparisons . 


(4.3) 


For  the  customer  removal  algorithm  A^  is  defined  by  relation  (2.31)  . 
By  changing  the  order  of  the  summation  in  (2.31),  we  can  conceptualize 

as  being  evaluated  while  the  k-th  busy  interval  progresses  and  requiring 
+  0(1)  storage  locations.  The  number  of  operations  required  are 

\(Mk  "  11 

- - -  comparisons , 

and  (4.4) 

M  (M  -  1) 

— — -  +  +  0(1)  additions. 


1 


In  the  time-contraction  algorithm  A^  is  defined  by  relation  (2.42). 
Since  A^  can  be  evaluated  as  the  busy  period  progresses  by  simply  adding 
service  requirements,  waiting  times,  and  arrival  times  relative  to  the 
start  of  the  busy  interval,  we  only  need  0(1)  storage  locations  and 
3M^  +  0(1)  additions  for  its  computation. 

From  the  above  discussion  we  see  that  of  all  three  algorithms,  the 
time-contraction  procedure  requires  no  buffer  with  randomly  varying  size. 

In  addition,  the  c  nputational  load  of  the  time-contraction  algorithm  is  ■ 
the  least  of  the  three.  Hence,  the  time-contraction  procedure  is  the  least 
costly  to  implement. 


5.  SIMULATION  RESULTS  AND  CONCLUSIONS 

Earlier  we  analyzed  the  asymptotic  bias  of  the  three  estimation 
algorithms.  Since  we  are  unable  to  calculate  the  variance  of  the  esti¬ 
mators  as  a  function  of  N,  the  question  of  consistency,  whether  the 
estimates  converge  asymptotically  to  the  delay  gradient,  remains  unan¬ 
swered.  In  addition,  for  an  M/D/1  queue  we  would  like  to  know  if  our 
algorithms  are  asymptotically  efficient,  whether  they  achieve  the  Cramer- 
Rao  bounds  derived  in  (3.28),  (3.57),  and  (3.69).  We  would  also  like 
to  investigate  the  robustness  of  the  customer-removal  and  time-contraction 
estimation  procedures  by  seeing  how  they  perform  for  a  variety  of  queueing 
systems.  We  attempt  to  answer  these  questions  by  presenting  the  results 
of  simulating  all  three  algorithms  for  an  M/D/1  queue  and  simulating  the 
customer-removal  and  time-contraction  algorithms  for  M/M/1,  D/M/1,  and 


U/M/l  queues. 
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We  simulate  a  single-server  queue  by  the  following  recursion  for 
successive  waiting  times: 


w  , ,  =■  max (w  +  x  -  0  /0}  with  w,  =  0  , 
n+1  n  n  n  J  1 


(5.1) 


where  xn  and  9^  are  random  variables  corresponding  to  the  n-th  service 
requirement  and  the  inter-arrival  time  between  the  n-th  and  (n+l)-st 
customer  respectively.  When  w  goes  to  zero,  this  signals  the  start 
of  a  new  busy  period. 

We  now  describe  the  calculations  necessary  to  evaluate  the  statistics 
of  a  given  estimator.  To  derive  estimates  for  the  mean  and  variance  <->f  the 

A 

delay  derivative  estimate  based  on  k  busy  periods,  D'  .  ,  we  generate  a 


certain  sample  size  of  k-busy  period  records ,  processing  each  to  form 

a  delay  gradient  estimate  D'  .  .  We  compute  estimates  of  the  bias 

(k) 

b(k)  and  variance  a ^  associated  with  6'  as  follows: 


i  s 

6  =■  —  y  6’ 

(£)  N  l  (k)  ,i 
s  i*-1 


3D 

3A  ' 


(.5.2) 


The  measure  of  performance  we  are  most  interested  in  is  the  fractional 


r.m.s.  error: 


/ec^-d'  )2 

l31  (kr 


C5.4) 
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We  approximate  f (k)  by  £(k)  defined  as 


£  (k) 


/^2  2 

v  b  (k)  +  tyOcn 

“  SD 

3X 


(5.5) 


The  value  for  Ng  employed  was  arrived  at  by  experimentation  as  Ng  =  400. 
For  each  type  of  queue  of  interest  we  computed  £ (k)  for  k  =  10,  100,  and 
1000  with  the  utilization  factor  -P  =  ^-x  -  varying  between  .1  and  .7. 

Having  defined  the  measure  f(k)  which  we  will  use  to  compare  the 
algorithms,  we  first  discuss  the  simulation  results  for  an  M/D/1  queue. 
Curves  of  f(k)  for  k  =  100  and  1000  are  presented  in  graphs  (5.1)  and 
(5.2) .  The  consistency  of  the  three  procedures  is  suggested  by  the  im- 

A 

proved  performance  with  increasing  N.  Together  with  f (k) ,  we  display 
the  fractional  r.m.s.  errors  that  follow  from  our  Cramer-Rao  bounds 
for  the  variance  of  each  algorithm.  The  closeness  of  the  simulation 
curves  to  their  respective  lower  bounds  suggest  that  all  three  algorithms 
are  asymptotically  efficient  for  an  M/D/1  queue. 

The  bounds  displayed  in  graphs  (5.1)  and  (5.2)  are  reflections  of  the 
amount  of  information  about  the  delay  gradient  contained  in  the  observa¬ 
tions  which  each  algorithm  employs.  In  actual  practice,  the  rate  of 
message  arrivals  to  the  queues  associated  with  each  link  in  a  Computer- 
Network  would  be  perturbed  periodically  by  some  global  routing  procedure. 
Assuming  that  the  adjustments  in  arrival  rates  are  small  and  that  we 
have  available  an  old  estimate  for  the  delay  gradient,  the  number  of  busy 
periods  of  observations  required  to  achieve  a  given  r.m.s.  estimation 
error  may  be  significantly  less  than  that  predicted  by  graphs  (5.1)— (5.2) . 


Fractional  RMS  error 


82933AW008 


Graph  (5.1)  Lower  Bounds  on  Fractional  RMS  Error  for 

m|d|  l(k-100) . 


Fractional  RMS  error 
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Utilization  factor-^ 

Graph  (5.2)  Lower  Bounds  on  Fractional  PMS  Error  forM|D|l. 
(k=*1000)  . 


Now  we  examine  the  robustness  properties  of  the  customer  removal 
and  time-contraction  algorithms  by  comparing  their  performance  for  M/M/1, 
U/M/l*  and  D/M/1  queues.  The  simulation  results  for  the  fractional  r.m.s. 
error  of  the  two  algorithms  for  M/M/1,  D/M/1,  and  D/M/1  queues  are  presented 
in  graphs  (5. 3) -(5. 8),  respectively.  The  time  contraction  algorithm 
performs  slightly  better  than  the  customer-removal  algorithm  for  M/M/1 
and  both  perform  nearly  the  same  for  a  U/M/l  queue.  The  only  dramatic 
differences  occur  for  the  D/M/1  queue,  where  the  time-contraction  pro¬ 
cedure  performs  better  than  the  customer-removal  algorithm.  This  result 
is  reasonable  since  for  a  D/M/1  queue  the  time-contraction  algorithm 
simulates  a  change  in  arrival  rate  in  the  exact  way  dictated  by  the 
structure  of  the  queue. 

Hence,  we  have  three  estimation  algorithms  that  appear  from  simulation 
and  analytical  results  to  be  asymptocially  unbiased,  consistent,  and  effi¬ 
cient  in  the  case  of  an  M/D/1  queue.  In  evaluating  the  most  promising 
algorithm  as  far  as  robustness,  computational  complexity,  and  storage 
requirements ,  we  can  only  choose  between  the  customer-removal  and  time- 
contraction  algorithms ,  since  the  customer-addition  procedure  as  formulated 
is  only  applicable  to  queues  where  all  customers  have  identical  service 
requirements.  The  simulation  results  suggest  the  time-contraction  algo¬ 
rithm  as  the  best  candidate  since  it  appears  to  be  at  least  as  robust 
as  the  customer-removal  procedure,  while  having  considerably  less  computa¬ 
tional  and  storage  requirements. 
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Graph  (5.3)  Time-Contraction  Algorithm  -  Fractional  RMS 
Errors  for  m|m|1. 
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Graph  (5.4)  Customer- Remcrval  Algorithm  -  Fractional  RMS 
Errors  for  m|m |l. 
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Graph (5.5)  Time-Contaction  Algorithm-  Fractional  RMS 
Error  for  u|m|1. 
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Graph (5.6)  Cus tome Removal  Algorithm  -  Fractional  RMS 
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Graph (5.7)  Time-Contraction  Algorithm  -  Fractional  SMS  Error 
for  D | M | 1 . 
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Graph (5.8)  Customer- Removal  Algorithm  -  Fractional  RMS 
Error  for  D | M | 1 . 
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